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Abstract 

We have used a novel computer with a novel integra- 
tion technique to study the evolution of the entire planetary 
system for nearly 100 million years. This calculation confirms 
that the evolution of the Solar System as a whole is chaotic, 
with a remarkably short timescale of exponential divergence 
of about 4 million years. Additional numerical experiments 
indicate that the dynamical evolution of the Jovian planets 
is chaotic apart from the rest of the system, as is the motion 
of Pluto. 

Advances in computer technology have made it possible to begin to di- 
rectly address the age-old question of the nature of the long-term evolution of 
the Solar System, with startling results. Sussman and Wisdom [1] presented 
numerical evidence that the motion of Pluto is chaotic, with a timescale for 
exponential divergence of nearby trajectories of only about 20 million years. 
Subsequently, Laskar [2] found numerical evidence of the chaotic evolution 
of the whole Solar System excluding Pluto, with a timescale for exponential 
divergence of only about 5 million years. Laskar 's calculation was feasible be- 
cause he analytically averaged the equations of motion to remove the rapid 
variations with timescales of order the orbital period. The averaged equations 
are perturbative and necessarily truncated after a particular order in eccentric- 
ity, inclination and mass ratio. A new integration of the whole Solar System 
without these approximations was clearly required. 

Direct integrations of the whole planetary system are computationally ex- 
pensive. Notable long-term integrations of the outer Solar System include: 
the classic 1 million year integration of Cohen, Hubbard, and Oesterwinter [3], 
the 5 million year integration of Kino shit a and Nakai [4], the 210 million 
year integration performed on the Digital Orrery [5], the 100 million year 
integration of the LONGSTOP project [6], and the 845 million year Digi- 
tal Orrery integration of Sussman and Wisdom [1]. Studies of the long-term 
evolution of the whole Solar System have been more limited because the com- 
putational resources required are significantly larger, by about two orders of 
magnitude. Integrations of the whole Solar System include: the 3 million year 
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Digital Orrery integration (excluding Mercury) [5], the 2 million year integra- 
tion of Richardson and Walker [7], and the recent ±3 million year integration 
of Quinn, Tremaine, and Duncan [8, 9] (hereafter QTD). 

We have developed new computational techniques and new computer hard- 
ware to make possible a direct integration of the whole Solar System spanning 
a significantly longer interval than previously achieved. Our new direct inte- 
gration of the equations of motion spans 36,000,000,000 days, or about 98.6 
million years. Our earlier result concerning the chaotic motion of Pluto, as well 
as the result of Laskar that the Solar System is chaotic are both confirmed. In 
order to localize the sources of the chaotic behavior we have carried out nu- 
merous additional long-term integrations. We have found that the evolution 
of the Jovian planets is independently chaotic, as is the motion of Pluto. 

Method of Integration 

We use the symplectic n-body mapping method of Wisdom and Holman [10] 
to integrate the planetary system. This new mapping method is a generaliza- 
tion of the mapping method of Wisdom [11, 12]. The basic idea is as follows. 

The Hamiltonian for the planetary n-body problem can be written 

H = H Kepler + H Interactions (1) 

where the first term represents the Keplerian motion of each of the planets with 
respect to the sun, and the second term describes the planetary perturbations. 
The simplest form of the map is obtained by adding short period terms so that 
the Hamiltonian becomes 

Hltfap = HKepler + H Interaction^ ^2w(^t)s (2) 

where S 2w (t) is a periodic sequence of Dirac delta functions with period 27T, 
and ft is the mapping frequency. The averaging principle [13] suggests that the 
additional short period terms do not significantly affect the long-period evolu- 
tion. This Hamiltonian is locally integrable: between the delta functions the 
motion is purely Keplerian and integration across the delta functions is easily 
carried out in Cartesian coordinates. The resulting map of the phase space 
onto itself is a composition of canonical transformations and is thus canonical 
(i.e., the Jacobian is a symplectic matrix). In general, to make a mapping of 
this sort each part of the Hamiltonian, considered by itself, must be not only 
integrable but also efficiently solvable. The Keplerian phase of the evolution 
can be efficiently carried out directly in Cartesian coordinates with the aid of 
Gauss's / and g functions. The uniform use of the Cartesian coordinates also 
eliminates any need for costly intermediate coordinate transformations. 

The mapping method can be made accurate to arbitrarily high-order in the 
mapping step-size. This and other refinements are discussed in Wisdom and 
Holman [10]. In the applications presented here the second order version of the 
n-body mapping is used. Second order is achieved by evolving the system with 
the Kepler Hamiltonian for a half mapping step, followed by an alternating 



succession of full interaction kicks and whole Keplerian steps, but ending with 
a half Kepler step. 

This ridiculously simple scheme is actually a remarkably good and efficient 
integrator. Wisdom and Holman [10] used this new map to compute the 
evolution of the outer planets for a billion years and compared the results to 
the results of the 845 million year integrations performed on the Digital Orrery 
by Sussman and Wisdom [1] using conventional integration techniques. All of 
the results of that study were confirmed, including details of the very long- 
period variations (> 500 million years) in the inclination of Pluto, and the 
20 million year exponential divergence of trajectories, which indicated chaotic 
behavior in the motion of Pluto. The new map is about an order of magnitude 
faster than traditional methods of integration. 

Our 100 Million Year Integration 

Our 100 million year integration of the planetary system was performed us- 
ing the Supercomputer Toolkit [14]. The Supercomputer Toolkit is the succes- 
sor to the Digital Orrery [15]. It is a small multiprocessor computer optimized 
for the numerical solution of systems of ordinary differential equations. The 
Toolkit was built as a collaboration between MIT and Hewlett-Packard. It is a 
family of standard software and hardware modules that can be interconnected 
in a variety of configurations as appropriate for particular applications. Each 
processor of the Toolkit is three times faster than the entire Digital Orrery, 
as measured by running the same computation. We used the eight-processor 
configuration of the Toolkit at MIT, with the new symplectic n-body map, to 
carry out eight 100 million year integrations of the planetary system. Each 
processor was used to run a separate integration with slightly different initial 
conditions so that we could look for exponential divergence. 

The most complete of the long-term integrations of the whole Solar Sys- 
tem is that of QTD. Particular care was taken by QTD to make their physical 
model accurate. They included general relativistic corrections, and a carefully 
crafted quadrupole approximation to account for the long-term effects due to 
the finite size of the Earth-Moon system. They used initial conditions and 
masses derived from JPL ephemeris DE102. Our physical model is the same 
as that of QTD except in our treatment of the effects of general relativity. 
General relativistic corrections can be written in Hamiltonian form, but we 
have not been able to integrate them analytically. Instead we used the po- 
tential approximation of Nobili and Roxburgh [16], which is easily integrated, 
but only approximates the relativistic corrections to the secular evolution of 
the shape and orientation of the orbit. 

Multistep methods typically require more than a hundred steps per orbit for 
stability; the mapping method is stable with as few as 10 steps per orbit. The 
step-size for our integration was rather arbitrarily chosen to be 7.2 days; with 
7.2-day steps output points are easily compared to the ephemeris of QTD. We 
integrated backward in time. A 22 million year integration using 3.6-day steps 
was carried out to check that the 7.2-day integration was sufficiently accurate. 



All of the position calculations were done in pseudo-quadruple precision, in 
an effort to reduce the effect of roundoff error. We believe now that the use 
of quadruple precision was an unnecessary precaution, and needlessly doubled 
the computation time of our experiment. Even on this demanding task each 
processor evolved the Solar System at the rate of about 30 years per second. 
Thus our eight 100 million year integrations took a total of about 1000 hours 
of Toolkit time. 

We have compared our integration to the reversed-time segment of the 
0.5-day step-size integration of QTD, and the agreement is quite good. The 
maximum difference in the argument of perihelion of Mercury over the 3 mil- 
lion year interval is of order 0.0001 radians; for comparison, the precession of 
the argument of perihelion due to general relativity is about 2w radians over 
3 million years. Evidently, our more approximate treatment of relativity is of 
little consequence. As another comparison, Table 1 lists the maximum differ- 
ences between the eccentricities of the planets compared to QTD. Also listed 
are the differences between the integration of QTD and that of Laskar [9]. The 





|Laskar - QTD| | 


Toolkit - Q' 


Mercury 


0.0041 


0.000018 


Venus 


0.0020 


0.000065 


Earth 


0.0024 


0.000059 


Mars 


0.0041 


0.000132 


Jupiter 


0.0038 


0.000047 


Saturn 


0.0081 


0.000162 


Uranus 


0.0051 


0.000008 


Neptune 


0.0026 


0.000002 


Pluto 


— 


0.000001 



Table 1: The maximum differences in the eccentricities of the planets in the 
integrations of Laskar, QTD, and the Toolkit show excellent agreement among 
the integrations. Laskar did not include Pluto in his calculation. 



table illustrates that the evolution computed with the mapping agrees quite 
well with the more conventional direct integration of QTD. The Toolkit inte- 
gration and QTD are mutually more consistent than either is to the integration 
of Laskar, though it is not clear whether this discrepancy is due primarily to 
model differences or to the approximations used by Laskar. The model dif- 
ferences are actually to our benefit because they help address the important 
question of the sensitivity of our results to model parameters. 

The Evolution of the Solar System is Chaotic 

Exponential divergence of nearby trajectories is indicative of chaotic behav- 
ior. The divergence of trajectories may have both quasiperiodic and exponen- 
tial components. Of course, if present, the exponential component eventually 



dominates. In order to estimate the Lyapunov exponent, we extracted the ex- 
ponential part of the divergence from the quasiperiodic part by the following 
procedure. We first converted the positions and momenta into Keplerian ele- 
ments, then formed the standard shape and orientation variables h = esintu, 
k = e cos tZ7, p — sin t/2 sin fi, and q = sin i/2 cos H, where e is the eccentricity, 
w is the longitude of perihelion, i is the inclination, and H is the longitude 
of the ascending node. The full set of these variables for all of the planets 
constitutes the "full secular phase space." Distance in this space is the or- 
dinary Euclidean distance. The quasiperiodic component is not as strong in 
the secular phase space as it is in the full phase space, allowing us to study 
the exponential part of the divergence more easily. If exponential divergence is 
observed in the secular phase space exponential divergence will also be present 
in the full phase space. 

The divergence in the secular phase space between two of the 100 million 
year calculations is shown in Figure 1. The initial conditions for the two calcu- 
lations differed by about 1 millimeter in the x coordinate of Pluto. The secular 
divergence has two distinct segments. In the latter segment the divergence is 
dominated by an exponential with a timescale of about 4 million years. The 
initial segment is apparently dominated by a smaller exponent with a timescale 
of about 12 million years. 

If the system is chaotic then the divergence of individual planets will also 
be exponential with the same exponent as the whole system if the trajectories 
are followed for sufficiently long time. Over shorter intervals the divergence 
of individual planets can be different, and examination of the individual di- 
vergences can give insight into the mechanisms responsible for the chaotic 
behavior in the system. 

In the inner Solar System the individual planet divergences in the secu- 
lar phase space are similar to that of the full secular phase space, with two 
distinct segments. In the outer Solar System, over most of the 100 million 
years spanned by our integrations the divergences appear to be dominated 
by a 12 million year exponential component, with evidence of the 4 million 
year component appearing in only the last 5 million years. A change in the 
timescale of exponential divergence could occur if the system went from one 
region of the phase space characterized by one exponent into another region 
characterized by the other exponent. This possibility can be ruled out because 
the 4 million year component appears much later in the outer planets than in 
the inner planets. A viable interpretation is that there are two distinct mech- 
anisms generating exponential divergence that are simultaneously operating. 
Apparently the inner planets are more sensitive indicators of the 4 million year 
process than are the outer planets. 

The evolution of Pluto over the full 100 million years has characteristics 
quite similar to those found in long term integrations of the outer planets. 
For example, we observe the 34 million year amplitude modulation of the 3.8 
million year oscillation of the argument of perihelion of Pluto [5] . 

Figure 2 shows the exponential divergence of the difference of positions of 
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Figure 1: Exponential divergence of nearby trajectories is indicated by the 
average linear growth of the logarithm of the full secular phase space distance 
between two different runs with very slight differences in initial conditions. The 
segment following the initial transient has an exponential timescale of about 
12 million years. The divergence is subsequently dominated by an exponential 
with a timescale near 4 million years. 
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Figure 2: The divergence of the distance between the positions of Pluto in two 
different runs with very slight differences in initial conditions is characterized 
by a timescale near 12 million years. 



Pluto in two of our 100 million year runs. These runs differed by 1 part in 10 15 
in one coordinate of the initial position of Mars. The plot is consistent with 
an exponential divergence of about 12 million years. This rate of exponential 
divergence is similar to that observed in the outer planets, and to the slower 
component observed in the full secular phase space. One interpretation is 
that Pluto is passively driven by chaotic behavior of the rest of the Solar 
System. Another interpretation is that Pluto is independently chaotic with a 
rate of exponential divergence which is only coincidentally similar to that of 
the rest of the system. We present evidence for the latter interpretation below. 
Considering the typical slow convergence of Lyapunov estimates, a 12 million 
year timescale of exponential divergence is in satisfactory agreement with the 
inverse Lyapunov exponent of 15 to 20 million years for Pluto reported by 
Sussman and Wisdom [1]. 

Secular Resonances 

Laskar [17, 18] has found three resonance angles which alternately librate 
and circulate. This may be an important corroboration of the chaos as indi- 
cated by the exponential divergence. The Laskar angles are 

^K-t^-ffti-n;), (3) 

<r 2 = 2K-«7;)-(n;-n;), (4) 

<Ts = K-w;)-(nj-n;). (5) 

The quantities zu° and fi° are the angles of the proper modes as denned by 
Laskar [17]. On the basis of the alternate libration of the angles a 2 and 03, 
which cannot simultaneously librate, Laskar concluded that the mechanism 
responsible for the chaotic behavior of the Solar System was resonance overlap 
of the corresponding secular resonances. 

In our calculation, <J\ and 02 also alternately circulate and librate, but 03 
just circulates, as shown in Figure 3. Our angles track Laskar 's angles for the 
initial segment of the computations, but they soon diverge enough so that the 
changes between libration and circulation occur at very different times. Such 
differences are consistent with the chaotic character of the evolution. The 
higher-order angle 

€T 4 = 3(«7j - «r;) - 2(n; - n|) (6) 

is also incompatible with 02 and 03 . In our calculation 04 has intervals of 
libration. The data presented by Laskar suggest that this is also the case in 
his calculation, though the angle is not explicitly mentioned. Since 1/1 < 
3/2 < 2/1 it is possible that we and Laskar are seeing different portions of 
a chaotic zone that spans the phase space from the 1/1 resonance to the 2/1 
resonance. In addition to these angles we also find that the angle 

<r 5 = (mi - w») + (n; - at) (7) 

alternates between circulation and libration. 
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Figure 3: In our calculation two of the angles of Laskar, a\ and <r 2 , alternately 
circulate and librate. Laskar*s third angle, 03, circulates. The higher order 
angle cr 4> which is incompatible with both cr 2 and 03, has intervals of li brat ion. 
The angle (7 5 also alternately circulates and librates. 



Though our calculations are completely consistent with those of Laskar, we 
are not fully convinced that his proposed mechanism accounts for the observed 
chaotic behavior of the Solar System. First, it is not clear that the alternate 
libration and circulation of (r 2 is indicative of a dynamically significant chaotic 
separatrix. In Figure 4 we show a polar plot of the amplitude and phase of 
the combination of proper modes corresponding to a\ and a 2 . Our polar plot 
of (T\ is similar to that shown in Laskar [17]. However, our polar plot of <r 2 
differs from that shown in Laskar [17] and Laskar, Quinn, and Tremaine [9]. 
For the polar amplitude they use just the amplitude of eccentricity mode 4, 
which is never close to zero. The amplitude corresponding to a 2 has additional 
factors, some of which approach zero. Using just the amplitude of eccentricity 
mode 4 gives an artificial impression of a transition from libration to circula- 
tion. Our plot, which includes all the factors in the amplitude, does not give 
a clear indication of a chaotic separatrix. Rather, it gives the impression of 
a complex high-dimensional trajectory projected onto a plane (see Figure 4). 
The alternation of libration and circulation of the polar angle may just be 
an artifact of the origin being in the midst of this complex projected trajec- 
tory. A simple integrable model for motion near the 2/1 commensurability in 
the planar-elliptic restricted three-body problem [19, 20] illustrates how the 
alternate circulation and libration of an angle can be misleading. There the 
resonance angle corresponding to the largest term in the disturbing function 
can show alternate circulation and libration, even though there is no chaotic 
behavior. The phenomenon can be eliminated by an appropriate change of 
variables. The polar plot corresponding to a\ is more like that expected of 
an angle associated with a chaotic mechanism. It not only alternately circu- 
lates and librates, but as it circulates it loiters around an apparent unstable 
equilibrium. There is also an excluded region near the center of the plot. 

Furthermore, there are too many unrelated angles which alternately cir- 
culate and librate. It might be expected that only one incompatible set of 
angles would show alternate circulation and libration. However, we see the 
phenomenon in angles involving unrelated modes. Eccentricity and inclina- 
tion modes 3 and 4 are involved in 02, 03, and 04; eccentricity modes 1,5, and 
8, and inclination modes 1, 2, and 8 are involved in <r\ and <t 5 . The two sets 
of modes are disjoint, and yet there are correlations in the behavior of angles 
associated with unrelated sets of modes. Most striking is the transition in be- 
havior in four of the angles near an integration time of 50 million years. This 
suggests that a single mechanism is driving all of the angles. If the mechanism 
is associated with one of the angles presented, we feel the most convincing is 
<T\. It is also possible that the mechanism generating the chaos is unrelated 
to all of these angles, and that they are all just sensitive indicators of chaotic 
irregularity of the underlying system trajectory. 

We have not identified any other angles whose motion is suggestive of a 
dynamical mechanism. The LONGSTOP team speculated [21] that the angle 
2w\ — 2zn? + Hy — ft| might alternately circulate and librate, and thereby 
provide evidence of chaotic behavior of the outer planets. We find that this 
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Figure 4: The variables are X{ = A{ cos <T{ and yi = Ai sin <7j, where .A,- are the 
amplitudes of the combination of the Laskar proper modes corresponding to 
the angles (7». The plot corresponding to a?, gives the impression of a complex 
high- dimensional trajectory projected onto a plane. The plot corresponding 
to <T\ shows some indication of a dynamically significant chaotic separatrix. 
There is a suggestion of an unstable equilibrium and associated asymptotic 
trajectories, and there is a region near the center that is avoided. 
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angle rotates uniformly with a period somewhat greater than 50 million years. 
This angle is one of 30 angles presented by the LONGSTOP team in their 
"numerology table." We find all 30 angles rotate uniformly; there is no hint 
of chaotic irregularity in the behavior of any of the angles. 

That particular secular resonances are responsible for the chaotic behavior 
of the Solar System could be established by an analytical demonstration that 
the resonances involved are sufficiently strong and close to give resonance 
overlap. At the moment, however, we feel no dynamical mechanism for the 
observed chaotic behavior of the Solar System has been clearly demonstrated. 

Chaotic Evolution of the Jovian Planets 

Despite the growing number of long-term integrations of the outer plan- 
ets, to date no integration has directly tested whether the massive planets 
themselves evolve chaotically. In our 845 million year integration of the outer 
planets [1], the orbital elements of Neptune appeared to have discrete line spec- 
tra, in marked contrast to the clearly broad-band character of the spectrum 
of Pluto. On the basis of this spectral evidence we dismissed the possibility 
of chaotic behavior in the Jovian planets. To be thorough, we have now car- 
ried out a new billion year evolution of the outer planets, using the mapping 
method, with a slightly perturbed initial position of Neptune. We found, to our 
surprise, that the subsequent evolution of the Jovian planets diverged expo- 
nentially from the first calculation, with a timescale of exponential divergence 
of only 5 million years! 

Our initial reaction was that there must be something wrong with our 
method of integration. To check this we carried out two new integrations 
of the outer planets spanning 100 million years using the traditional linear 
multistep Stormer predictor, the same integrator used in the Digital Orrery 
integrations. The initial conditions and masses were the same as in the Dig- 
ital Orrery integrations. The integrations were carried out in ordinary IEEE 
double precision (64 bits) with a step-size of 32.7 days, the same step-size 
used in the 845 million year Digital Orrery integrations. We found that the 
trajectories of the Jovian planets diverged exponentially with a timescale of 
about 19 million years. Apparently, we were misled by the spectral evidence. 

To further check that this result did not depend on either the step-size 
or the precision of the calculation we carried out four more integrations of 
the outer planets using the Stormer predictor. Each spanned more than 400 
million years. In these integrations the accumulation of position was carried 
out in pseudo-quadruple precision, as in the Digital Orrery integrations. One 
pair used a step-size of 32.7 days, the special Digital Orrery step-size; the other 
pair used an arbitrarily chosen step-size of 28 days. The initial conditions were 
the same as in our earlier outer planet integrations. In one run of each pair 
the initial position of Neptune was perturbed by 7.5mm. The energy errors 
again grew linearly with time, with slopes between that of our 210 million year 
Digital Orrery integration and that of our 845 million year integration. Both 
pairs of runs gave remarkably consistent results. The secular divergence of the 
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Jupiters in the 28-day runs is shown in Figure 5. The Jovian planets diverged 
exponentially with a timescale somewhat longer than 20 million years. 

The map and Stormer calculations both indicate that the motion of the 
Jovian planets is chaotic. However, they are discrepant in the estimate of 
the exponential timescale. We believe this is a result of the fact that the 
trajectory computed by the mapping method is not exactly the same as the 
trajectory determined by the initial conditions. The mapping differs from the 
actual n-body dynamics by the addition of extra high-frequency terms. The 
averaging principle suggests that these high-frequency terms do not contribute 
to the long-term evolution, and the close agreement of our results with those of 
Laskar once again confirms the validity of the use of averaging. However, the 
initial conditions used in an averaged system should properly take into account 
the presence of the extra high-frequency components in the unaveraged system. 
The use of the same initial conditions in the averaged and unaveraged systems 
corresponds to slightly different initial conditions for the long- term evolutions. 
With the mapping we do not yet know how to properly adjust the initial 
conditions to account for the extra high-frequency terms. In our outer-planet 
integrations these slight uncontrolled adjustments of effective initial conditions 
appear to yield different estimates of the Lyapunov exponent. 

To investigate this further we carried out 8 integrations of the outer plan- 
ets using the map with different step-sizes, ranging from 0.3 years to 1 year. 
Changing the step-size changes the high-frequency components, and slightly 
changes the effective initial condition for the long-term evolution. Each inte- 
gration spanned about 300 million years. We found that the measured diver- 
gence timescale varied from about 3 million years to about 30 million years. 
The dispersion in the estimates of the Lyapunov exponent are much larger 
than the dispersion observed in the Stormer runs. The Lyapunov exponent 
was not obviously correlated with step-size, in particular the estimate of the 
Lyapunov exponent was not monotonic with step-size. In one of these runs, 
with step-size near 0.617979 years, the motion of the outer planets was clearly 
quasiperiodic; the secular phase space divergence did not grow exponentially. 
The results suggest that the Lyapunov exponent for the Jovian planets is not 
a simple function of the initial conditions. Most nearby initial conditions lead 
to exponential divergence (most with a shorter timescale for exponential di- 
vergence than that obtained with the Stormer integrations), but there are also 
nearby initial conditions that do not give chaotic behavior. 

With each outer planet integration we ran a pair of massless Plutos, with 
initial positions differing by about 1cm. A remarkable result is that the expo- 
nential divergence of Plutos always has a timescale between 10 and 20 million 
years, independent of how chaotically the Jovian planets behave. This is true 
even in the most extreme runs, where the Jovian planets were quasiperiodic, 
and where the Jovian planets diverged with a timescale of 3 million years. This 
clearly demonstrates that the mechanism generating the chaotic behavior in 
the motion of Pluto is extremely robust, and independent of whether the rest 
of the system is chaotic. 
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Figure 5: The secular phase space divergence of the trajectories of Jupiter in 
the 28 day Stormer integrations show a timescale for exponential divergence 
that is somewhat longer than 20 million years. The divergence saturates after 
about 250 million years at a rather small value, perhaps indicating a hidden 
constraint on the trajectories. 
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Numerical Chaos? 

The fact that almost all long-term integrations of the Solar System give 
exponential divergence of trajectories with a timescale in the range of 3-30 
million years in physically quite different models is very striking and unsettling. 
Another unsettling feature of the chaotic behavior we observe in long-term 
planetary integrations is that nothing dramatic happens. This is compounded 
by the fact that in no case have we unambiguously identified the mechanism 
responsible for the chaotic behavior. This lack of mechanism and lack of 
obvious irregular behavior is in marked contrast to the clearly understood 
mechanisms and irregular character observed in other examples of chaotic 
behavior in the solar system, for example, asteroids on chaotic trajectories 
near commensurabilities [11, 12, 22], the chaotic tumbling of Hyperion [23] 
and other irregularly shaped satellites [24], and the chaotic motion of the 
Uranian satellites near commensurabilities [25, 26, 27]. 

Perhaps the exponential divergence is a numerical artifact? The detailed 
agreement of the billion year evolution of Pluto using different integrators is 
impressive evidence against this perverse possibility. Furthermore, the detailed 
agreement between our 100 million year Solar System integration and that of 
Laskar is particularly convincing, because of the radically different methods 
used. 

To further convince ourselves that not all long-term integrations are subject 
to some universal numerical instability, which yields a spurious exponential di- 
vergence, we have carried out a control integration. We have integrated the 
outer planets without Uranus for about 250 million years, with the mapping 
method. Over that period we see no sign of exponential divergence of the re- 
maining Jovian planets. This integration, together with the isolated quasiperi- 
odic integration mentioned above, shows that numerical models of planetary 
systems can, in fact, show quasiperiodic behavior on a several hundred million 
year time-scale. Long-term integrations do not always give positive Lyapunov 
exponents. 

Altogether, the evidence for the chaotic behavior in these long-term plan- 
etary integrations is very convincing, but there remains the logical possibility 
that the exponential divergence is a subtle numerical artifact. To positively 
conclude that the chaos observed in these long-term planetary integrations 
is not a result of numerical artifacts requires an unambiguous identification 
of a physical mechanism and an analytic evaluation to determine that the 
mechanism actually accounts for the observed chaos. 

Conclusions and Speculations 

Our 100 million year integration of the entire Solar System indicates that 
the Solar System is chaotic with a timescale for exponential divergence of 
about 4 million years. The fact that we find similar behavior in all respects 
to the calculation of Laskar strongly supports the conclusion that the Solar 
System is chaotic. That we and Laskar have carried out different kinds of 
numerical experiments, with slightly different masses, slightly different initial 
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conditions, and slightly different physics, shows that the chaotic character of 
the Solar System is not sensitively dependent on the precise model or numerical 
methods. 

Our experiments indicate that the Jovian planets by themselves behave 
chaotically for most initial conditions near our reference system, though our 
estimates of the Lyapunov exponent have rather large dispersion. 

All of our estimates of the Lyapunov exponent of Pluto give approximately 
the same divergence timescale of 10-20 million years, with different methods 
of integration, different planetary masses, different initial conditions, and even 
independently of whether the rest of the system is behaving chaotically. Our 
earlier result that the evolution of Pluto is chaotic is thus multiply confirmed. 

We will not fully understand the consequences of the observed chaotic evo- 
lution of the Solar System until we clearly understand the dynamical mecha- 
nisms responsible for it. Though in our calculation the behavior of the secular 
resonance angles is consistent with those in Laskar's calculation, the identifi- 
cation of resonance overlap of particular secular resonances as the mechanism 
generating the chaotic behavior of the Solar System is not unambiguously 
demonstrated. Our numerical experiments suggest that there are at least two 
independent mechanisms generating chaotic behavior. One mechanism gen- 
erates chaos in the Jovian subsystem, and Pluto is independently chaotic in 
the field of the Jovian planets. Yet another mechanism is suggested by the 
simultaneous presence of two different exponential timescales in our full So- 
lar System integrations. Secular resonances among the inner planets may be 
driving the faster timescale, as Laskar suggested. However, the most convinc- 
ing of the secular resonances involve Mercury. Both Mercury and Pluto have 
high eccentricity and inclination, which strongly couples the eccentricity and 
inclination secular subsystems. Perhaps one of the mechanisms generating the 
chaos originates with Mercury, and is similar to the mechanism generating the 
chaos in the motion of Pluto? 

References 

[1] G.J. Sussman and J. Wisdom Science 241, 433 (1988). 

[2] J. Laskar Nature 338, 237 (1989). 

[3] C.J. Cohen, E.C. Hubbard, and C. Oesterwinter, in Astronomical Papers 
of the American Ephemeris (Government Printing Office, Washington, 
DC, 1973), vol. 22, pp. 1-92. 

[4] H. Kinoshita, and H. Nakai Celestial Mechanics 34, 203 (1984). 

[5] J. Applegate, M.R. Douglas, Y. Giirsel, G.J. Sussman, and J. Wisdom 
Astron. J., 92, 176 (1986). 



16 



[6] A.E. Roy, I.W. Walker, A.J. Macdonald, LP. Williams, K. Fox, CD. 
Murray, A. Milani, A.M. Nobili, P.J. Message, A.T. Sinclair, and M. 
Carpino, Vistas in Astronomy 32, 95 (1988). 

[7] D.L. Richardson, and C.F. Walker, in Astrodynamics 1987 (Advances 
in the Astronautical Sciences, Vol. 65), J.K. Soldner, A.K. Misra, R.E. 
Lindberg, and W. Williamson, eds. (Univelt, San Diego, 1987), 1473. 

[8] T.R. Quinn, S.D. Tremaine, and M. Duncan Astron. J. 101, 2287 (1991). 

[9] J. Laskar, T.R. Quinn, and S.D. Tremaine, Icarus 95, 148 (1992). 

[10] J. Wisdom, and M. Holman Astron. J. 102, 1528 (1991). 

[11] J. Wisdom Astron. J. 87, 577 (1982). 

[12] J. Wisdom Icarus 56, 51 (1983). 

[13] V.I. Arnold Mathematical Methods of Classical Mechanics (Springer- 
Verlag, 1978). 

[14] H. Abelson, A.A. Berlin, J. Katzenelson, W. McAllister, G.J. Rozas, G.J. 
Sussman, and J. Wisdom, preprint. 

[15] J. Applegate, M.R. Douglas, Y. Gursel, P. Hunter, C. Seitz, and G.J. 
Sussman, IEEE Trans. Comput. C34(1985). 

[16] A. Nobili and I. Roxburgh, in Relativity in Celestial Mechanics and As- 
trometry, Kovalevsky, Brumberg eds. (Reidel, Dordrecht, 1986). 

[17] J. Laskar Icarus 88, 266 (1990). 

[18] J. Laskar in Proc. IAU Symposium 152, Ferraz-Mello ed. (1991). 

[19] J. Wisdom Celestial Mechanics 38, 175 (1986). 

[20] J. Henrard, A. Lemaitre, A. Milani, and CD. Murray Celestial Mechanics 
38, 335 (1986). 

[21] A. Nobili, A. Milani, and M. Carpino Astron. Astrophys. 210, 313 (1989). 

[22] J. Wisdom Icarus, 63, 272 (1985). 

[23] J. Wisdom, S.J. Peale, and F. Mignard Icarus 58, 137 (1984). 

[24] J. Wisdom Astron. J. 94, 1350. 

[25] W. Tittemore and J. Wisdom Icarus 74, 172 (1988). 

[26] W. Tittemore and J. Wisdom Icarus 78, 63 (1989). 

[27] W. Tittemore and J. Wisdom Icarus 85, 394 (1990). 

17 



[28] This work depended upon many people. The Supercomputer Toolkit com- 
puter was designed and built with the help and generous support of the 
Hewlett Packard Company. We especially thank Joel Birnbaum of Hewlett 
Packard for his personal support. William McAllister led the team at HP. 
Dan Zuras of HP helped with the development of the scientific subrou- 
tines we used. At MIT Andy Berlin, Bill Rozas, Henry Wu, Hal Abelson, 
and Jacob Katzenelson contributed to the hardware and software de- 
sign. Matthew Holman helped with the development and testing of the 
mapping integrator. We thank Scott Tremaine, Jacques Laskar, Thomas 
Quinn, and Martin Duncan for helpful discussions. 

This report describes research done at the Artificial Intelligence Labora- 
tory of the Massachusetts Institute of Technology. Support for the labora- 
tory's artificial intelligence research is provided in part by the Advanced 
Research Projects Agency of the Department of Defense under Office of 
Naval Research contract N00014-89-J-3202 and in part by the National 
Science Foundation under grant number MIP-9001651. Jack Wisdom is 
also supported in part by a NSF Presidential Young Investigator Award 
AST-887365 and in part by the NASA Planetary Geology and Geophysics 
Program under grant NAGW-706. We are grateful for the hospitality of 
the California Institute of Technology where both authors are on leave 
from MIT. 



18 



